% 我看见电磁场于空间穿行...二次元电磁场的时域有限差分（FDTD）模拟
% 可能有bug
% Gitee Repo


clc
clear

dx = 0.1;
dt = 0.01;
TICK = 1000;
L=20;
[x y] = meshgrid(0:dx:L);
N = size(x,1);

Ez_2 = zeros(N, N);
Bx_2 = zeros(N+1, N);
By_2 = zeros(N, N+1);
Ez_1 = zeros(N, N);
Bx_1 = zeros(N+1, N);
By_1 = zeros(N, N+1);

jz = 500;

figure()

function j = current(t)
		j = 1/(exp(-5*t+5)+1);
		if j > 0.99
				j = 1;
		end
    %	T = 10;
    %	theta = 2*pi/T*t;
    %	if theta < pi/2
    %			j = sin(theta);
    %	else
    %			j = 1;
    %	end
end

for tick = 1:TICK
  % Yee grid
  Bx_2(2:N,:) = Bx_1(2:N,:) - dt/dx*(Ez_1(2:N,:) - Ez_1(1:N-1,:));
  By_2(:, 2:N) = By_1(:,2:N) + dt/dx*(Ez_1(:,2:N) - Ez_1(:,1:N-1));
  Ez_2(:,:) = Ez_1(:,:) - dt/dx*(Bx_2(2:N+1,:) - Bx_2(1:N,:)) + dt/dx*(By_2(:,2:N+1) - By_2(:,1:N));

  % Source: Current
	Ez_2(uint32(N/2), uint32(N/2)) -= dt*jz*current(tick*dt);

  % Source: alternating electric field E
  % Ez_2(uint32(N/2),uint32(N/2)) = 2*sin(2*pi/15*tick*dt);

  % Boundary
  Bx_2(1,:)= 0;
  Bx_2(N+1,:)= 0;
  Bx_2(:,1)= 0;
  Bx_2(:,N)= 0;

  By_2(1,:)= 0;
  By_2(N,:)= 0;
  By_2(:,1)= 0;
  By_2(:,N+1)= 0;

  % Update
  Bx_1 = Bx_2;
  By_1 = By_2;
  Ez_1 = Ez_2;

  if mod(tick,20) == 0
    % Spatial Average (Since B is defined on half grid)
    _Bx_2 = (Bx_2(2:N+1,:) + Bx_2(1:N,:))/2;
    _By_2 = (By_2(:, 2:N+1) + By_2(:,1:N))/2;
    _Bx_1 = (Bx_1(2:N+1,:) + Bx_1(1:N,:))/2;
    _By_1 = (By_1(:, 2:N+1) + By_1(:,1:N))/2;
    % Time Average (Since B is defined on half timestep)
    __Bx = (_Bx_2 + _Bx_1)/2;
    __By = (_By_2 + _By_1)/2;

    __Ez = Ez_1;

    __zeros = zeros(N);

    clf
    hold on
    axis equal
    axis([0 L 0 L -2 2])
    view([30 45])
%    view([0 90])
    caxis([-1 1])

		wire = line([x(uint32(N/2), uint32(N/2)) x(uint32(N/2), uint32(N/2))], [y(uint32(N/2), uint32(N/2)), y(uint32(N/2), uint32(N/2))], [-5 5],'color','r');
		set(wire,'color',[current(tick*dt) 0 0]);

    DRAW_RES = 10;
    _range = 1:DRAW_RES:N;
    quiver3(x(_range,_range), y(_range,_range), __zeros(_range,_range), __zeros(_range,_range), __zeros(_range,_range), __Ez(_range,_range), 0);
    quiver3(x(_range,_range), y(_range,_range), __zeros(_range,_range), __Bx(_range,_range), __By(_range,_range), __zeros(_range,_range), 0);
%    surf(x,y,(__Ez.^2 + __Bx.^2 + __By.^2),'Edgecolor','none');
%    any(any(__Ez.^2 - (__Bx.^2 + __By.^2) > 0.1))
    drawnow
    pause(0.1)
%    print(['EM2D_PDE_' num2str(tick) '.png'])

  end

end
